با استفاده از داده های زلزله ها در ایران و جهان به سوالات زیر پاسخ دهید.
بارگزاری داده ها و کتابخانه ها:
library(readr)
library(dplyr)
library(ggmap)
library(highcharter)
library(stringr)
earth_q <- read_rds("data/historical_web_data_26112015.rds")
۱. با استفاده از داده های historical_web_data_26112015.rds و استفاده از نمودار پراکنش سه بعدی بسته plotly نمودار طول، عرض و عمق زلزله ها را رسم نمایید. علاوه بر آن بزرگی هر نقطه را برابر بزرگی زمین لرزه قرار دهید.
library(plotly)
p <- plot_ly(data = earth_q, x = ~Latitude, y = ~Longitude, z = ~Depth, size = ~Magnitude,
marker = list(color = ~Magnitude, symbol = 'circle', sizemode = 'diameter', colorscale = c('#FFE1A1', '#683531'), showscale = TRUE), sizes = c(0.5, 7.8),
text = ~paste('Province:', Province, '<br>City:', City, '<br>Magnitude:', Magnitude)) %>%
add_markers() %>%
layout(scene = list(title = 'Iran Earthquakes between Sep. to Nov. 2015',
xaxis = list(title = 'Latitude', range = c(23, 40)),
yaxis = list(title = 'Longitude', range = c(40, 71)),
zaxis = list(title = 'Depth', range = c(0, 162))))
p۲. پویانمایی سونامی های تاریخی را بر حسب شدت بر روی نقشه زمین رسم نمایید.(از داده زلزله های بزرگ استفاده نمایید.)
برای این سوال از داده های بزرگ استفاده می کنیم، سپس مقادیر نامعلوم را داده حذف کرده و بر اساس مقادیر شدت، نمودار را رسم می کنیم.
library(gganimate)
disaster = read_delim("data/disaster.txt", "\t", escape_double = FALSE, trim_ws = TRUE) %>%
rename(lat = LATITUDE,long = LONGITUDE, magnit = INTENSITY,name = COUNTRY,year = YEAR) %>%
dplyr::select(lat, long, magnit, name, year)
nadis <- na.omit(disaster) %>% arrange(-magnit)
mapWorld <- borders("world", colour="gray50", fill="gray50") # create a layer of borders
p <- ggplot() + xlab("longitute") + ylab("latitude") + mapWorld + ggtitle("Earthquake Intensity")
p <- p + geom_point(aes(x = nadis$long, y = nadis$lat, size = nadis$magnit, frame = nadis$magnit), color="indianred1") + guides(size=guide_legend(title="Intensity"))
gganimate(p)
animation <- gganimate(p, "images/eq.gif")۳. نمودار چگالی دو بعدی زلزله های تاریخی ایران را رسم کنید.( از داده iran_earthquake.rds و لایه stat_density_2d استفاده نمایید).
برای حل این سوال ابتدا یک داده ی پرت را که خارج از ایران بود حذف می کنیم. سپس نمودار چگالی را بر روی نقاط رخداد زلزله می کشیم.
iran_eq <- read_rds("data/iran_earthquake.rds")%>% arrange(-Long)
iran_eq = iran_eq[-c(1),]
p <- ggplot(iran_eq, aes(x=Long, y=Lat)) +
geom_point() + stat_density_2d(aes(fill = ..level..), geom = "polygon") +
xlab("Longitude") + ylab("Latitude") +
guides(fill=guide_legend(title="Density")) +
scale_fill_distiller(palette=4, direction=-1) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
ggtitle("Iran Earthquakes Density Plot")
p۴. احتمال اینکه در ایران در پنج سال آینده زلزله به بزرگی هفت ریشتر رخ دهد را محاسبه کنید. (از احتمال شرطی استفاده کنید.)
iran_sig_eq <- iran_eq %>% filter(Mag >= 7) %>% arrange(desc(OriginTime))
last <- iran_sig_eq %>% slice(1)
diff = as.numeric(Sys.time() -last$OriginTime, units="days")
diff = diff %/% 365
diff_year = diff + 5
after_year = 1 # for first one
for (i in 1:nrow(iran_sig_eq)) {
if (i > 1) {
differ = as.numeric(iran_sig_eq[i-1,]$OriginTime - iran_sig_eq[i,]$OriginTime, units="days") %/% 365
if (differ >= diff_year) {
after_year = after_year + 1
}
}
}
after_year = after_year / nrow(iran_sig_eq)
after_two = 1 # for first one
for (i in 1:nrow(iran_sig_eq)) {
if (i > 1) {
differ = as.numeric(iran_sig_eq[i-1,]$OriginTime - iran_sig_eq[i,]$OriginTime, units="days") %/% 365
if (differ >= 2) {
after_two = after_two + 1
}
}
}
after_two = after_two / nrow(iran_sig_eq)
p = (after_year)/(after_two)
cat("Probability of earthquake in 5 years is", p)Probability of earthquake in 5 years is 0.5
۵. بر اساس داده های زلزله های بزرگ ابتدا تعداد و متوسط کشته زلزله ها را بر حسب کشور استخراج نمایید. سپس نمودار گرمایی تعداد کشته ها را بر روی کره زمین رسم نمایید.(مانند مثال زیر!)
برای حل این سوال ابتدا داده های NA را حذف می کنیم، سپس میانگین و مجموع تلفات زلزله ها را بدست می آوریم. برای نمایش کشورها بر روی نقشه نیاز است که نام کشورها را به نام کوتاه استاندارد آن ها تبدیل کنیم. به همین دلیل از کتابخانه ی countrycode استفاده می کنیم. سپس نمودار گرمایی را یک بار برای تعداد کشته ها و یک بار برای میانگین تعداد کشته ها توسط کتابخانه ی rworldmap می کشیم.
library(rworldmap)
library(countrycode)
disaster_sum <- read_delim("data/disaster.txt", "\t", escape_double = FALSE, trim_ws = TRUE) %>%
rename(lat = LATITUDE,long = LONGITUDE, magnit = INTENSITY,country = COUNTRY,year = YEAR, death=TOTAL_DEATHS) %>%
dplyr::select(lat, long, magnit, country, year, death)
disaster_sum = na.omit(disaster_sum)
disaster_sum <- disaster_sum %>% group_by(country) %>% summarise(fatality_sum = sum(death), fatality_mean = mean(death))
disaster_sum$country <- countrycode(disaster_sum$country, "country.name", "iso3c", warn = TRUE, nomatch = NA,
custom_dict = NULL, custom_match = NULL, origin_regex = FALSE)
knitr::kable(disaster_sum)| country | fatality_sum | fatality_mean |
|---|---|---|
| AFG | 934 | 155.66667 |
| ALB | 2735 | 390.71429 |
| DZA | 19021 | 1118.88235 |
| ARG | 516 | 103.20000 |
| ARM | 187000 | 62333.33333 |
| AUS | 12 | 12.00000 |
| AUT | 1 | 1.00000 |
| AZE | 310087 | 77521.75000 |
| PRT | 1122 | 561.00000 |
| BGD | 200 | 200.00000 |
| BEL | 2 | 2.00000 |
| BIH | 1 | 1.00000 |
| BGR | 138 | 34.50000 |
| CAN | 28 | 28.00000 |
| CHL | 70059 | 1893.48649 |
| CHN | 1964956 | 11101.44633 |
| COL | 2420 | 242.00000 |
| CRI | 99 | 24.75000 |
| HRV | 5005 | 1668.33333 |
| CYP | 42 | 21.00000 |
| DJI | 2 | 2.00000 |
| DOM | 8 | 4.00000 |
| ECU | 111202 | 22240.40000 |
| EGY | 1138 | 379.33333 |
| SLV | 1328 | 265.60000 |
| ETH | 70 | 35.00000 |
| FRA | 650 | 650.00000 |
| GEO | 280 | 93.33333 |
| GHA | 25 | 12.50000 |
| GRC | 29047 | 660.15909 |
| GLP | 5006 | 2503.00000 |
| GTM | 23054 | 4610.80000 |
| GIN | 443 | 443.00000 |
| HND | 3 | 3.00000 |
| ISL | 11 | 11.00000 |
| IND | 58776 | 5877.60000 |
| IDN | 4481 | 128.02857 |
| IRN | 637473 | 20563.64516 |
| IRQ | 70200 | 17550.00000 |
| ISR | 35100 | 11700.00000 |
| ITA | 353786 | 7075.72000 |
| JPN | 6151 | 307.55000 |
| KAZ | 451 | 225.50000 |
| KGZ | 79 | 39.50000 |
| MKD | 1083 | 361.00000 |
| MEX | 10997 | 845.92308 |
| MNG | 30 | 30.00000 |
| MNE | 131 | 131.00000 |
| MAR | 16728 | 5576.00000 |
| MMR | 503 | 251.50000 |
| NPL | 11771 | 3923.66667 |
| NLD | 1 | 1.00000 |
| NZL | 259 | 64.75000 |
| PAK | 146342 | 24390.33333 |
| PNG | 687 | 98.14286 |
| PER | 84300 | 2218.42105 |
| PHL | 13425 | 383.57143 |
| PRT | 82319 | 11759.85714 |
| ROU | 18 | 6.00000 |
| RUS | 12025 | 1717.85714 |
| LCA | 900 | 900.00000 |
| SRB | 6 | 2.00000 |
| SVN | 22 | 11.00000 |
| SLB | 62 | 20.66667 |
| ZAF | 1 | 1.00000 |
| ESP | 2907 | 969.00000 |
| CHE | 9 | 9.00000 |
| SYR | 105700 | 35233.33333 |
| TWN | 10035 | 501.75000 |
| TJK | 15529 | 5176.33333 |
| THA | 1 | 1.00000 |
| TTO | 1 | 1.00000 |
| TUN | 48000 | 24000.00000 |
| TUR | 806913 | 19680.80488 |
| TKM | 110001 | 55000.50000 |
| UKR | 11 | 11.00000 |
| USA | 902 | 27.33333 |
| USA | 168 | 84.00000 |
| UZB | 4880 | 4880.00000 |
| VEN | 43808 | 5476.00000 |
| YEM | 4060 | 1015.00000 |
heatMap <- joinCountryData2Map(disaster_sum, joinCode = "ISO3",
nameJoinColumn = "country")80 codes from your data successfully matched countries in the map
1 codes from your data failed to match with a country code in the map
165 codes from the map weren't represented in your data
mapCountryData(heatMap, nameColumnToPlot = "fatality_sum",
mapTitle="Earthquakes total fatality", oceanCol=gray(0.3),
colourPalette=c("cadetblue1","cadetblue3", "deepskyblue2", "deepskyblue3", "deepskyblue4", "dodgerblue4", "navy"), missingCountryCol = "white")heatMap <- joinCountryData2Map(disaster_sum, joinCode = "ISO3",
nameJoinColumn = "country")80 codes from your data successfully matched countries in the map
1 codes from your data failed to match with a country code in the map
165 codes from the map weren't represented in your data
mapCountryData(heatMap, nameColumnToPlot = "fatality_mean",
mapTitle="Earthquakes fatality rate", oceanCol=gray(0.3),
colourPalette=c("cadetblue1","cadetblue3", "deepskyblue2", "deepskyblue3", "deepskyblue4", "dodgerblue4", "navy"), missingCountryCol = "white")۶. با استفاده از داده لرزه های بزرگ و به وسیله طول، عرض، شدت، عمق مدلی برای پیش بینی تعداد کشته های زلزله بیابید.
برای مدلسازی از glm استفاده می کنیم. به علت آنکه توزیع ما توزیع پسین حاشیه ای دارد که نشان می دهد واریانس توزیع نرمال نامعلوم است، از خانواده ی inverse gamma استفاده می کنیم.
disaster_death <- read_delim("data/disaster.txt", "\t", escape_double = FALSE, trim_ws = TRUE) %>%
rename(lat = LATITUDE,long = LONGITUDE, magnit = INTENSITY, depth = FOCAL_DEPTH, country = COUNTRY,year = YEAR, death=TOTAL_DEATHS) %>%
dplyr::select(lat, long, magnit, depth, country, year, death)
# generalized regression model
glm_model <- glm(death ~ lat + long + magnit + depth, family = Gamma(link = "inverse"), data = disaster_death)
summary(glm_model)
Call:
glm(formula = death ~ lat + long + magnit + depth, family = Gamma(link = "inverse"),
data = disaster_death)
Deviance Residuals:
Min 1Q Median 3Q Max
-4.055 -3.273 -2.611 -1.587 7.681
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.158e-03 2.561e-04 4.524 7.90e-06 ***
lat -9.875e-07 2.824e-06 -0.350 0.7267
long -2.387e-07 4.959e-07 -0.481 0.6304
magnit -1.070e-04 2.406e-05 -4.445 1.12e-05 ***
depth 7.511e-06 3.933e-06 1.910 0.0568 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Gamma family taken to be 14.9034)
Null deviance: 3965.6 on 429 degrees of freedom
Residual deviance: 3398.8 on 425 degrees of freedom
(5596 observations deleted due to missingness)
AIC: 5834
Number of Fisher Scoring iterations: 10
glm_model <- glm(death ~ magnit + depth, family = Gamma(link = "inverse"), data = disaster_death)
summary(glm_model)
Call:
glm(formula = death ~ magnit + depth, family = Gamma(link = "inverse"),
data = disaster_death)
Deviance Residuals:
Min 1Q Median 3Q Max
-3.895 -3.270 -2.619 -1.607 7.898
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.223e-03 2.326e-04 5.259 2.30e-07 ***
magnit -1.153e-04 2.197e-05 -5.247 2.44e-07 ***
depth 6.468e-06 1.967e-06 3.287 0.0011 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Gamma family taken to be 15.12215)
Null deviance: 3965.6 on 429 degrees of freedom
Residual deviance: 3409.4 on 427 degrees of freedom
(5596 observations deleted due to missingness)
AIC: 5832.1
Number of Fisher Scoring iterations: 10
cat("Test GLM model using null and model deviances: ",1-pchisq(3965.6 - 3409.4, df=(429 - 427)))Test GLM model using null and model deviances: 0
مشاهده می کنیم که متغیر های طول و عرض جغرافیای تاثیری ندارند و عمق و شدت زلزله تاثیر زیادی بر روی تعداد کشته ها دارد. هم چنین با توجه به خروجی مشاهده می کنیم که خطای ما حول ۲ است. از طرفی آخرین خروجی به ما نشان می دهد که مدل ما از مدل intercept-only بهتر عملکرد و فرض صفر که یکسان بودن این مدل است باطل است.
۷. با استفاده از داده worldwide.csv به چند سوال زیر پاسخ دهید. تحقیق کنید آیا می توان از پیش لرزه، زلزله اصلی را پیش بینی کرد؟
worldwide <- read_csv("data/worldwide.csv")۸. گزاره " آیا شدت زلزله به عمق آن بستگی دارد" را تحقیق کنید؟ (طبیعتا از آزمون فرض باید استفاده کنید.)
برای این سوال از تست کوریلیشن و تست استقلال chi استفاده می کنیم.
cor.test(worldwide$depth, worldwide$mag, method = 'spearman')
Spearman's rank correlation rho
data: worldwide$depth and worldwide$mag
S = 2.8439e+13, p-value < 2.2e-16
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
0.2344252
eq_matrix <- worldwide %>% select(depth, mag) %>% filter(depth > 0) %>% as.matrix()
chisq.test(eq_matrix)
Pearson's Chi-squared test
data: eq_matrix
X-squared = 534440, df = 59889, p-value < 2.2e-16
در تست کوریلیشن مشاهده می کنیم که احتمال اینکه این دو متغیر از دو توزیع غیرمرتبط آمده باشند، صفر است. پس فرض صفر باطل است و این دو متغیر به یکدیگر وابسته هستند. همچنین در تست استقلال نیز مشاهده می کنیم که احتمال استقلال دو متغیر بسیار کم بوده و در نتیجه فرض صفر باطل است و این دو متغیر از یکدیگر مستقل نیستند. پس فرض عدم ارتباط عمق و شدت زلزله رد می شود.
۹. میانگین سالانه زلزله ها را بر حسب کشور به دست آورید. آیا میتوان دلیلی در تایید یا رد تئوری هارپ ارائه کرد.
برای حل این سوال، ابتدا می بایست کشور محل رخداد زلزله را بدست بیاوریم. برای این کار از نقاط نقشه ی rworldmap استفاده می کنیم. سپس بر اسا کشور و سال گروهبندی کرده و تعداد زلزله ها و میانگین شدت زلزله ها را بدست می آوریم. در نهایت نمودار موزاییک زلزله ها را نمایش می دهیم.
library(sp)
coords2country = function(points)
{
countriesSP <- getMap(resolution='high')
#setting CRS directly to that from rworldmap
pointsSP = SpatialPoints(points, proj4string=CRS(proj4string(countriesSP)))
# use 'over' to get indices of the Polygons object containing each point
indices = over(pointsSP, countriesSP)
indices$ISO3 # returns the ISO3 code
}
s <- worldwide %>% select(longitude, latitude)
s$country <- coords2country(s)
s <- s %>% select(country)
worldwide <- bind_cols(worldwide, s)
worldwide <- worldwide %>% filter(!is.na(country)) %>%
mutate(year = as.numeric(format(time, format = "%Y")))
worldwide_sum <- worldwide %>% group_by(year, country) %>%
summarise(count = n(), mean_mag = mean(mag), mean_depth = mean(depth))
hchart(worldwide_sum %>% filter(year == 2015), "treemap", hcaes(x = country, value = count, color = mean_mag)) %>%
hc_title(text = "2015 Earthquakes", style = list(fontWeight = "bold"))hchart(worldwide_sum %>% filter(year == 2016), "treemap", hcaes(x = country, value = count, color = mean_mag)) %>%
hc_title(text = "2016 Earthquakes", style = list(fontWeight = "bold"))hchart(worldwide_sum %>% filter(year == 2017), "treemap", hcaes(x = country, value = count, color = mean_mag)) %>%
hc_title(text = "2017 Earthquakes", style = list(fontWeight = "bold"))hchart(worldwide_sum %>% filter(year == 2018), "treemap", hcaes(x = country, value = count, color = mean_mag)) %>%
hc_title(text = "2018 Earthquakes", style = list(fontWeight = "bold"))پروژه ی هارپ از سال ۲۰۱۵ منحل شده است. در صورتی که این پروژه تاثیر زیادی بر روی زلزله داشته باشد، باید با کاهش چشمگیر زلزله مواجه باشیم. اما مشاهده می کنیم که زلزله ها افزایش داشته اند. همچنین روند واضحی در این میان دیده نمی شود. در نتیجه نمی توانیم در مورد تاثیر پروژه ی هارپ نظری دهیم.
۱۰. سه حقیقت جالب در مورد زلزله بیابید.